function md = madist(data)
sd = size(data, 1);
m = (1 / sd) .* sum(data, 1);
data = data - m(ones(sd, 1), :);
[Q,R] = qr(data,0);

ri = R' \ data';
md = sum(ri .* ri, 1)' * (sd - 1);